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ABSTRACT 

We derive the first secondary eclipse map of an exoplanet, HD 189733b, based on Spitzer IRAC 
8 micron data. We develop two complementary techniques for deriving the two dimensional planet 
intensity: regularized slice mapping and spherical harmonic mapping. Both techniques give similar 
derived intensity maps for the infrared day-side flux of the planet, while the spherical harmonic method 
can be extended to include phase variation data which better constrain the map. The longitudinal 
offset of the day-side hot spot is consistent with that found in prior studies, strengthening the claim 
of super-rotating winds, and eliminating the possibility of phase variations being caused by stellar 
variability. The latitude of the hot-spot is within 12.5° (68% confidence) of the planet's equator, 
confirming the predictions of general circulation models for hot Jupiters and indicative of a small 
planet obliquity. 

Subject headings: Methods: observational — Techniques: photometric — Infrared: planetary systems 
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1. INTRODUCTION 

One of the great challenges in studying extrasolar plan- 
ets is that we cannot directly resolve the surfaces of these 
bodies. This is a common problem in astronomy, and one 
solution is to use occultations or eclipses to spatially re- 
solve astronomical sources. Eclipse mapping has been 
appli e d variously to map stars (e.g. iSchmitt fc Kurster 
1993; iLanza et al.l Il998| ). accretion disks (e.g. Horne 
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1985f ). planetary sat ellites (e .g. Aksncs & Franklin 1975i), 
dwarf planets (e.g. lYoung e t al. 1999, 2001), and unre- 
solved radio sources (e.g. pTavlor T966, .1967) . Here we 
present the first eclipse map of an extrasolar planet. 

The secondary eclipse is the decrease in flux that occurs 
when a planet passes behind its host star. The eclipse 
ingress and egress contain information about the spatial 
distribution of the planet's day side fiux. For example, 
a planet with a brighter western hemisphere has an ap- 
parent delay in the times of ingress and egress, and a 
modified shape of the ingress and egress. This delay in 
eclipse timing due to a zonally asymmetrical planet was 
pred icted by \y illiam s ct al. (2006) and first observed by 
iKiTutson et all pOOTf ) and lAgol et al.l (|2010( ). In addition 
to this artificial time offset, the ingress/egress morphol- 
ogy may be us e d to map exoplanets as pointed out in 
iRauscher et"all (|2007f ). 

Longitudinal variations in the brightness of the planet 
cause a variation in the brightness of the planet with or- 
bital phase. T his orbital modulatio n has been observed 
in transiting (iKnutson et al.l [20071 l2 009b . a) an d non- 
trans iti ng systems (iCowan et al.l l2007i : iCrossfield et al.l 
I2010D . ICowan fc Agoll (|2008[ ) described how to invert 
such phase variations into a low-resolution longitudinal 
map of a planet; th e first such map was presented in 
IKnutson et all (|2007l ). 
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The short-period t ransiting jovian planet HD 189733b 
(|Bouchv et all 120051 ) is one of the first planets to 
have been detec ted in infrared emi ssion during sec- 
ondary eclipse ([Deming et al.l I2006D . and has been 
studi e d in depth since tha t detection (IKnutson et al.l 
20071 lGrillma.ir et all l2008t ICha rbonnc au et al.l 120081: 
Knutson et al.l l2009bl : lAgol et al.l l2010l: iDesert et al.l 

2qTW . 

We apply our mapping techn iques to th e seve n sec- 
ondary eclipse measurements of lAgol et al.l (|20100 and, 
in the case of spherical h armonic mapping, th e phase 
measurement presented in IKnutson et al.l (|2007l ). These 
data have been corrected for detector effects and stellar 
variability, and the seven eclipse light curves were binned 
in phase to 644 po ints. To this we appe nded the phase 
variation data from lKnutson et ahl (|2007[ ). binned by 275. 

2. ECLIPSE MAPPING 

2.1. Model Assumptions 

Three conditions need to be satisfied in order to con- 
struct an eclipse map: the planet emission pattern must 
be static, the time of superior conjunction must be 
known, and the planet limb-darkening must be negligi- 
ble. 

A static diurnal temperature pattern is required for 
any form of thermal mapping. Most general circulation 
models (GCMs) for hot Jupiters indicate < 1% vari- 
ability (Showman & Guillot 2002; Cooocr & ShowmM 
[20051 [Showman et al.l 120091 : IDobbs-Dixon et al.l \25im . 
while observations of multiple secondary eclipses of 
HD 189733b indicate that the eclipse depth does not vary 
significantly, < 2 7% over a period of more than 500 days 
(jAgol et al.ll2010( ). Hence the static map assumption is 
likely justified for this planet. 

Eclipse mapping can be applied provided that the time 
of superior conjunction is known to a few seconds, which 
is beyond the precision of radial velocity measureme- 
nents. We assume a circular orbit, hence pegging the 
time of eclipse to the time of transit which is precisely 
known. This assumption can then be checked by com- 
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paring the phase-function and secon dary eclipse map s, 
and seems to apply for HD 189733b (|Agol et al.ll201of ). 

Ignoring the effect of limb darkening is required so that 
the surface brightness of a location on the planet does 
not change as the planet rotates. This assumption is not 
required for our slice-mapping technique, since in that 
case we implicitly neglect changes in the planet's phase. 
Regardless, limb darkening is thought to be very weak for 
hot Jupiters, about 20% linear limb-darkening, especially 
in the wavelength of the observations (8 micron) , so this 
assumption is likely reasonable for HD 189733b, as we 
show below. 

In addition to these assumptions about the planet, we 
assume some things of the data. Any stellar variability 
must be corrected for, which may not be possible if it has 
a period close to that of the planet. We have corrected 
for stellar variability using the correlation betwee n opti- 
cal an d infrared flux variations as discussed in Agol et al.l 
(|2OT0t) . Detector ramp and oth er detector-related effects 
are corrected for as described in lAgol et al.l (|2Q1CI[ ). Addi- 
tionally, all light curves henceforth are taken to have the 
stellar flux subtracted off, so that measurements when 
the planet is completely occulted have zero value. 

We describe two approaches to planet mapping: (1) 
slice mapping (> j2.2[) : (2) spherical harmonic mapping 

2.2. The Slice Method 

During a time interval within ingress/egress, we can 
compute the area of the planet that disappears/appears 
behind the star if we assume a circular orbit. The frac- 
tion of flux that disappears during that time interval is 
then proportional to the mean surface brightness of the 
area that disappeared. This is the basis of our first map- 
ping technique. 

As a planet passes behind its star, the light curve 
flux, F(t), at any given time may be associated with the 
roughly crescent-shaped portion of its day side currently 
visible to the earth. It is thus proportional to the planet's 
light distribution function, I{6,(j)), integrated over the 
visible area, A{t). If we take two successive observa- 
tions at times tj-i and tj, the change in observed flux, 
AFj = F{tj) — F{tj-i), is due solely to the change in 
area AAj = A{tj) — A{tj^i) times the intensity averaged 
over that area. The average intensity, 7^ , over the slice 
with area AAj is then given by Ij = AFj / AAj . 

We partition the planet's ingress and egress into ob- 
servations at a finite number of intervals, Atj, and from 
these calculate each Ij, yielding two 'slice maps': ID 
maps of the planet day side. This process and its results 
for HD 189733b are shown in Figure [T] 

Once these maps were found, we combined them into 
a single 2-dimensional map of the planet, forming a grid 
of the ingress and egress slices. This step of the inver- 
sion requires a prior since if the ingress and egress maps 
are each made from N slices, the combined map contains 
roughly iV^ cells. We assume that the planet's intensity 
distribution varies smoothly, which is imposed via lin- 
ear regularization, in which a solution is chosen that is 
consistent with the data and has the smallest differences 
between adjacent slices._ We use a goodness of fit param- 
eter of -I- A J2ii^i ~ li-i)"^- We chose a regularization 
coefficient, A — 100, to give a reduced chi-square of or- 



der unity, resulting in the map of HD 189733b shown in 
Figure [11 

2.3. Spherical harmonic mapping 

Our second method for extracting a map from the 
eclipse light curve uses spherical harmonics. We applied 
this to two data sets: 1) only the secondary eclipse data; 
2) both the secondary eclipse and phase-function data. 

Taking 1(9, (j)) to be the intensity map of a transit- 
ing extrasolar planet, we may approximate it to an ar- 
bitrary degree of accuracy by truncating its harmonic 
sum at some finite number of terms, < I < Imax, 
with —l<m<l. We compute the light curve that 
would be observed if the planet's map consisted solely 
of each harmonic, and then solve for the coefficients of 
each harmonic light curve from the data via linear re- 
gression. We implemented the light curve calculation us- 
ing the HEALFix IDL package, which contains routines 
for the evaluation a nd projection of spherical harmonics 
(IGorski et al.l 120051) . The map is then formed by trans- 
forming back to the harmonics and summing the result. 
This transformation is portrayed graphically in Figure 
[21 and the results of its application to observations of 
HD 189733b can been seen in Figure [H This is essen- 
tially a two dimensional generalization of the Fourier de- 
composition te c hniqu e for phase mapping developed in 
Cowan fc A goT ("20081); note that the 2D image shown in 
Knutson et al. (200'^ only contained longitudinal infor- 
mation, while the latitudinal depdence was arbitrarily 
chosen. 

We generate the light curve for the real spherical har- 
monics, '^i (where i labels the harmonics used in our in- 
version) by painting a planet with a particular harmonic 
and then stepping through the planet's orbit for each of 
the observed orbital phases. At each phase the harmonic 
is integrated over the visible portion of the planet to get 
the total flux that would be observed. We integrate over 
the visible area by finding a function with curl equal to 
the harmonic and integrate it around the area's bound- 
ary, which is simply the concatenation of two known arcs, 
using Green's Theorem. This technique yields both the 
planet phase variation and secondary eclipse shape for 
each harmonic, as shown in Figure [51 

With a set of n spherical harmonic light curves, 
Ai(i), A„(t), we decompose the observed light curve, 
F{t), by computing the the coefficients ci, ...,c„ as Ci = 
J F{t)Xi{t)dt. These translate directly to the map's co- 
efficients for the harmonics V^i, 4'„. By adding the 
terms together we form a best-fit map of the entire 
planet, the closest approximation we can reach using the 
chosen harmonics, shown in Figure [4l 

For the HD 189733b Channel 4 IRAC map, we 
chose to use only the first four harmonics (/, m) = 
(0, 0), (1, — 1), (1, 0), (1, 1) which gave a reduced chi- 
square of the fit close to unity. Thus, the maps only 
contain a dipole variation in flux, which corresponds to 
the "hot spot" feature seen in simulations and inferred 
from the phase function. 

To quantify the statistical uncertainty on the measured 
hot spot location, we ran a Markov chain using the spher- 
ical harmonic fits. Beginning with a uniformly bright 
planet in a 4-dimensional vector space corresponding to 
the first four harmonic coefficients, we attempted jumps 
with magnitude based upon the uncertainty found in the 
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Fig. 1. — Top: Diagram demonstrating the slice mapping technique. Bottom: Ingress and egress maps (left and right), as well as combined 
map (center) of HD 189733b at 8 micron. The image is centered on the sub-observer point and the black regions are 50% as bright as the 
white. 

on the Metropolis-Hastings algorithm. The longitudi- 
nal and latitudinal positions of the hottest point on the 
planet model was recorded for each of the 25,000 itera- 
tions, and the resulting distribution is plotted in Figure 
H We found a hot spot latitude of 13.0 ± 14.1 deg North 
of the equator for the secondary eclipse inversion and 
3.1 ± 9.4 deg North for the eclipse plus phase function 
inversion. This indicates that the hot spot is consistent 
with being located at the equator. The corresponding 
longitude was 27.1^9^0^ deg East and 21.8 ± 1.5 deg East 
for the eclipse or eclipse plus phase inversion, respec- 
tively. Consequently, the phase function gives significant 
leverage in constraining both the longitude and latitude 
of the hot spot. However, we estimate that both are sub- 
ject to « 15% systematic errors due to neglect oi I > 1 
modes in the inversion based on our simulated data, as 
described next. 

3. APPLICATION TO SIMULATED DATA 

We created simulated light curves using a two dimen- 
sional map at c onstant pressure from a GCM for a hot 
Jupiter planet (jShowman et al.|[2008f ). to which we ap- 
plied the spherical harmonic mapping technique. We 
took a snap shot of the temperature at 100 mbar from 
a simulation, converted it to surface brightness in the 
Rayleigh- Jeans limit, and then created a simulated light 
curve from the images. The light curve was computed 
for identical orbital and size parameters as HD 189733a, b 
and for the same range of phase that was observed (about 
1/3 of the planet's orbit prior to secondary eclipse), and 
we assumed that the planet has a fixed temperature pat- 
tern with no limb-darkening, as well as a circular orbit 
and spherical bodies. We carried out six different ver- 
sions of the simulated light curves: 1) we first decom- 
posed the temperature map into spherical harmonics up 
to order Imax , and then used the map based on only these 
harmonics to compute a light curve with no noise; 2) we 
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Fig. 2. — A few spherical harmonics and their corresponding 
'harmonic light curves,' (Aq, Ai, A2, ...). 



initial decomposition for each coefficient, and a jump 
factor adjusted every 100 steps to achieve and accep- 
tance rate of « 44%, with the acceptance criterion based 




Degrees North of equator 

Fig. 3. — The probability distribution of the longitude and lat- 
itude of the location of the hot spot in the Imax = 1 spherical 
harmonic maps based on the secondary eclipse (solid lines) and 
secondary eclipse plus phase function (dashed lines). Angles are 
measured in the rotation frame of the planet, with the origin at 
the sub-stellar point. 

computed the same light curve, but with noise added a t 
the level of the observed data from lAgol et all (l20H : 
3) same as 2, but we shifted the hot spot North in lat- 
itude by 12, 23, 40 and 60 degrees; 4) we computed a 
light curve using the full temperature map; 5) we added 
noise to the light curve using the full map; 6) we carried 
out an identical test as version 1, but with 20% linear 
limb-darkening. We found that for Imax = 1, we could 
recover the correct coefficients of the spherical harmonics 
in version 1, which shows that the code performed as ex- 
pected. When noise was added, we found a scatter in the 
derived parameters at the 10% level. With Imax = 2, we 
found that the derived spherical harmonics became very 
poorly constrained, likely due to degeneracies in the har- 
monic light curves, so for the remainder of the analysis 
we used Imax — 1 when extracting the map. For version 
3, we ran 100 simulations at each latitudinal offset, and 
found that the mean of the recovered latitudes matched 
the input value, with a scatter varying from 12 deg RMS 
at the equator to 6 deg nearer to the pole. For version 4, 
we found that the amplitude of the spherical harmonics 
relative to the (0, 0) component had systematic errors as 
large as w 50%. However, the amplitude of the spheri- 
cal harmonics relative to the day-side flux inferred from 
the map agreed to better than 15%, while the location 
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Fig. 4. — A MoUeweide projection of the 2D infrared map of HD 
189733b with Imax = 1; latitude and longitude lines are spaced 
every 30°. Coordinates are centered on the sub-stellar point and 
intensity values are relative (white is brightest, black is 30.2% as 
bright) . 

of the hot spot differed by w 10 — 15% in longitude and 
latitude. For version 5 we found that the statistical error 
did not change much for Imax = 1, still at the 10% level. 
For version 5, we find that 20% limb darkening affects 
the derived coefficients by less than 5%. 

Based on these tests, we are confident in our inference 
of the longitude and latitude of the planet hot spot the 
dipole amplitude relative to the day-side flux at the w 
15% level. Additionally, we confirm the validity of our 
assumption to ignoring the effects of limb darkening for 
hot Jupiters in 8 micron. 

4. DISCUSSION AND CONCLUSIONS 
4.1. Slice vs. Spherical Harmonic Eclipse Maps 

The fact that the slice map (Figure [1]) and spherical 
harmonic map (Figure |4]) give qualitatively similar re- 
sults is reassuring as the two techniques have very differ- 
ent assumptions and potential systematic errors. Both 
maps favor a hot spot that is offset slightly latitudinally; 
however, when the phase function is included in the anal- 
ysis, this constrains the spot to be located more closely 
to the equator. This is due to the fact that the phase 
function provides a strong constraint on the longitudinal 
variation of the planet brightness, and thus does not al- 
low as much freedom in varying both the longitude and 
latitude to fit the secondary eclipse. 

The slice map has the disadvantage that it does not 
account for the rotation of the planet or the tilt of the 
planet, nor can it be used in conjunction with the phase 
variation of the planet, resulting in the north-south 
asymmetry visible in the slice map above absent from 
the harmonic mapping results. In addition, it is difficult 
to derive uncertainties on the intensity, which are highly 
correlated due to the regularization. 

The spherical harmonic method has the advantage of 
being able to easily account for these efi^ects, but the 
disadvantage that it assumes to a large degree that the 
planet intensity varies smoothly as we are constrained to 
use small Imax, in this case Imax = 1- 

4.2. Phase versus Eclipse Mapping 

Phase mapping (jCowan fc Agoll [20M ) offers more 
leverage on the longitudinal brightness distribution of 
the planet than does eclipse mapping. This is because 
the same change in flux is spread out over an orbital pe- 
riod, P, rather than the duration of ingress/egress, Ti„g, 
allowing for greater signal-to-noise constraints given the 
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same photometric uncertainties. The ratio of these two 
timescales is ri„g/P cx Rp/ina) («0.4% for HD 189733), 
where Rp and a are the planetary radius and orbital semi- 
major axis, respectively. Whenever it is feasible to obtain 
photometry throughout an entire orbit (or a sizable frac- 
tion thereof), phase curve inversion will provide the best 
estimate of the hottest and coldest local stellar times. 
This will be especially true as one tries to apply these 
techniques to smaller planets with longer periods. 

Eclipse mapping offers three significant advantages 
over phase mapping, however: 

1) Eclipse mapping is sensitive to to the latitudinal in- 
tensity distribution on a planet if it has non-zero impact 
parameter as it passes behind the disk of its star. 

2) The planet's brightness map is assumed to be static 
over the eclipse duration, TocI, which is much shorter than 
that assumed for phase mapping, P, and thus influenced 
less by other forms of variability. The ratio of these 
two timescales is Ted/P oc R*/a (« 3% for HD 189733). 
While there is no clear evidence yet for weather on tran- 
siting planets, it will certainly be the case that the in- 
tensity distributions of eccentric planets as seen by an 
observer will vary throughout an orbit. Eclipse mapping 
therefore offers the only prospect for constraining the in- 
tensity distributions of this intriguing class of objects. 

3) There is no theoretical limit to the spatial resolution 
one can achieve with eclipse mapping, although there 
may be some degeneracy as we have found. This is in 
stark contrast to phase mapping, where one runs into 
fundame ntal degeneracies wit h odd harmonics after only 
5 terms (jCowan fc Ago]|[2008l) . 

4.3. Summary 

We developed two techniques for converting an 
observed secondary eclipse light curve into a two- 
dimensional map of a planet. 

The spherical harmonic formalism is a robust way to 
make an exoplanet map from secondary eclipse light 
curves, and can accommodate observations of thermal 



phase variations. Slice maps may also be constructed 
at ingress and egress, then combined via regularization. 
While this technique seems to produce comparable day- 
side maps for the data at hand, it is not as robust or 
versatile as the spherical harmonic technique. 

We applied both mapping techniques to the transit- 
ing planet HD 189733b. We find that the primary 
day-side hotspot is offset to the East by 21.8(1.5)°, in 
agreement with the phase variation measurements of 
iKnutson et al.l ('20071 and the eclipse timing constraint 
from lAgol et al . (2010). The data indicate a hot spot 
latitude of 3.1(9.4)° North, consistent with an equato- 
rial hot spot. The fact that the hottest point lies close 
to the equator is consistent with a small obliquity of 
the planet as oblique rota tors are predicted to show 
strong latitudinal contrast (jLangton fc Laughlinl 12007): 
this conclusion is also consistent with theoretical compu- 
tations which predict a small obliquity for hot Jupiters 
(jFabrvckv et al.l 120071: iLevrard et all I2007D . Taken to- 
gether, these studies strongly favor super-rotating winds 
in the vicini ty of the 8 micron photosphere, as predicted 
in G CMs (CooDcr & Showman 2005; Showman et alj 

20ia 
The 



20081: [Showman et al 2009: R auscher fc Menou 
Burrows et at] l2010t IShowman fc Polvanil [20111 ). 



mapping codes used in this paper, written in IDL, are 
available from the authors upon request; however, the 
assumptions that go into this planet mapping technique 
need to be verified for other planets studied. 
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